/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
 *
 * 
 *                This source code is part of
 * 
 *                 G   R   O   M   A   C   S
 * 
 *          GROningen MAchine for Chemical Simulations
 * 
 *                        VERSION 3.2.0
 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
 * Copyright (c) 2001-2004, The GROMACS development team,
 * check out http://www.gromacs.org for more information.

 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 * 
 * If you want to redistribute modifications, please consider that
 * scientific software is very special. Version control is crucial -
 * bugs must be traceable. We will be happy to consider code for
 * inclusion in the official distribution, but derived work must not
 * be called official GROMACS. Details are found in the README & COPYING
 * files - if they are missing, get the official version at www.gromacs.org.
 * 
 * To help us fund GROMACS development, we humbly ask that you cite
 * the papers on the package - you can find them in the top README file.
 * 
 * For more info, check our website at http://www.gromacs.org
 * 
 * And Hey:
 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
 */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <math.h>
#include <sys/types.h>
#include <stdio.h>
#include <string.h>
#include <errno.h>
#include <ctype.h>
#include "futil.h"
#include "sysstuff.h"
#include "typedefs.h"
#include "smalloc.h"
#include "macros.h"
#include "gmxfio.h"
#include "txtdump.h"
#include "physics.h"
#include "macros.h"
#include "names.h"
#include "string2.h"
#include "symtab.h"
#include "gmx_fatal.h"
#include "warninp.h"
#include "vsite_parm.h"

#include "grompp.h"
#include "toputil.h"
#include "toppush.h"
#include "topdirs.h"
#include "gpp_nextnb.h"
#include "topio.h"
#include "topshake.h"
#include "gmxcpp.h"
#include "gpp_bond_atomtype.h"
#include "genborn.h"
#include "maths.h"

#define CPPMARK  	'#'	/* mark from cpp			*/
#define OPENDIR  	'['	/* starting sign for directive		*/
#define CLOSEDIR 	']'	/* ending sign for directive		*/

#define PM()  \
  printf("line: %d, maxavail: %d\n",__LINE__,maxavail()); \
  fflush(stdout)

static void free_nbparam(t_nbparam **param,int nr)
{
  int i;

  for(i=0; i<nr; i++)
    sfree(param[i]);
  sfree(param);
}

static int copy_nbparams(t_nbparam **param,int ftype,t_params *plist,int nr)
{
  int i,j,f;
  int nrfp,ncopy;

  nrfp = NRFP(ftype);
  
  ncopy = 0;
  for(i=0; i<nr; i++)
    for(j=0; j<=i; j++)
      if (param[i][j].bSet) {
	for(f=0; f<nrfp; f++) {
	  plist->param[nr*i+j].c[f] = param[i][j].c[f];
	  plist->param[nr*j+i].c[f] = param[i][j].c[f];
	}
	ncopy++;
      }
  
  return ncopy;
}

static void gen_pairs(t_params *nbs,t_params *pairs,real fudge, int comb, gmx_bool bVerbose)
{
    int     i,j,ntp,nrfp,nrfpA,nrfpB,nnn;
    real    scaling;
    ntp       = nbs->nr;
    nnn       = sqrt(ntp);
    nrfp      = NRFP(F_LJ);
    nrfpA     = interaction_function[F_LJ14].nrfpA;
    nrfpB     = interaction_function[F_LJ14].nrfpB;
    pairs->nr = ntp;
  
  if ((nrfp  != nrfpA) || (nrfpA != nrfpB))
    gmx_incons("Number of force parameters in gen_pairs wrong");

  fprintf(stderr,"Generating 1-4 interactions: fudge = %g\n",fudge);
  if (debug) {
    fprintf(debug,"Fudge factor for 1-4 interactions: %g\n",fudge);
    fprintf(debug,"Holy Cow! there are %d types\n",ntp);
  }
  snew(pairs->param,pairs->nr);
  for(i=0; (i<ntp); i++) {
    /* Copy param.a */
    pairs->param[i].a[0] = i / nnn;
    pairs->param[i].a[1] = i % nnn;
    /* Copy normal and FEP parameters and multiply by fudge factor */

    
    
    for(j=0; (j<nrfp); j++) {
        /* If we are using sigma/epsilon values, only the epsilon values 
         * should be scaled, but not sigma. 
         * The sigma values have even indices 0,2, etc.
         */
        if ((comb == eCOMB_ARITHMETIC || comb == eCOMB_GEOM_SIG_EPS) && (j%2==0))
            scaling = 1.0;
        else
            scaling = fudge;
        
        pairs->param[i].c[j]=scaling*nbs->param[i].c[j];
        pairs->param[i].c[nrfp+j]=scaling*nbs->param[i].c[j];
    }
  }
}

double check_mol(gmx_mtop_t *mtop,warninp_t wi)
{
  char    buf[256];
  int     i,mb,nmol,ri,pt;
  double  q;
  real    m;
  t_atoms *atoms;

  /* Check mass and charge */
  q=0.0;

  for(mb=0; mb<mtop->nmoltype; mb++) {
    atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
    nmol  = mtop->molblock[mb].nmol;
    for (i=0; (i<atoms->nr); i++) {
      q += nmol*atoms->atom[i].q;
      m  = atoms->atom[i].m;
      pt = atoms->atom[i].ptype;
      /* If the particle is an atom or a nucleus it must have a mass,
       * else, if it is a shell, a vsite or a bondshell it can have mass zero
       */
      if ((m <= 0.0) && ((pt == eptAtom) || (pt == eptNucleus))) {
	ri = atoms->atom[i].resind;
	sprintf(buf,"atom %s (Res %s-%d) has mass %g\n",
		*(atoms->atomname[i]),
		*(atoms->resinfo[ri].name),
		atoms->resinfo[ri].nr,
		m);
	warning_error(wi,buf);
      } else 
	if ((m!=0) && (pt == eptVSite)) {
	  ri = atoms->atom[i].resind;
	  sprintf(buf,"virtual site %s (Res %s-%d) has non-zero mass %g\n"
		  "     Check your topology.\n",
		  *(atoms->atomname[i]),
		  *(atoms->resinfo[ri].name),
		  atoms->resinfo[ri].nr,
		  m);
	  warning_error(wi,buf);
	  /* The following statements make LINCS break! */
	  /* atoms->atom[i].m=0; */
	}
    }
  }
  return q;
}

static void sum_q(t_atoms *atoms,int n,double *qt,double *qBt)
{
    double qmolA,qmolB;
    int     i;

    /* sum charge */
    qmolA = 0;
    qmolB = 0;
    for (i=0; i<atoms->nr; i++)
    {
        qmolA += atoms->atom[i].q;
        qmolB += atoms->atom[i].qB;
    }
    /* Unfortunately an absolute comparison,
     * but this avoids unnecessary warnings and gmx-users mails.
     */
    if (fabs(qmolA) >= 1e-6 || fabs(qmolB) >= 1e-6)
    {
        *qt  += n*qmolA;
        *qBt += n*qmolB;
    }
}

static void get_nbparm(char *nb_str,char *comb_str,int *nb,int *comb,
                       warninp_t wi)
{
  int  i;
  char warn_buf[STRLEN];
  
  *nb   = -1;
  for(i=1; (i<eNBF_NR); i++)
    if (gmx_strcasecmp(nb_str,enbf_names[i]) == 0)
      *nb = i;
  if (*nb == -1)
    *nb = strtol(nb_str,NULL,10);
  if ((*nb < 1) || (*nb >= eNBF_NR)) {
    sprintf(warn_buf,"Invalid nonbond function selector '%s' using %s",
	    nb_str,enbf_names[1]);
    warning_error(wi,warn_buf);
    *nb = 1;
  }
  *comb = -1;
  for(i=1; (i<eCOMB_NR); i++)
    if (gmx_strcasecmp(comb_str,ecomb_names[i]) == 0)
      *comb = i;
  if (*comb == -1)
    *comb = strtol(comb_str,NULL,10);
  if ((*comb < 1) || (*comb >= eCOMB_NR)) {
    sprintf(warn_buf,"Invalid combination rule selector '%s' using %s",
	    comb_str,ecomb_names[1]);
    warning_error(wi,warn_buf);
    *comb = 1;
  }
}

static char ** cpp_opts(const char *define,const char *include,
                        const char *infile,
                        warninp_t wi)
{
  int  n,len;
  int  ncppopts=0;
  const char *cppadds[2];
  char **cppopts  = NULL;
  const char *option[2] = { "-D","-I" };
  const char *nopt[2]   = { "define", "include" };
  const char *ptr;
  const char *rptr;
  char *buf;
  char warn_buf[STRLEN];
  
  cppadds[0] = define;
  cppadds[1] = include;
  for(n=0; (n<2); n++) {
    if (cppadds[n]) {
      ptr = cppadds[n];
      while(*ptr != '\0') {
	while((*ptr != '\0') && isspace(*ptr))
	  ptr++;
	rptr = ptr;
	while((*rptr != '\0') && !isspace(*rptr))
	  rptr++;
	len = (rptr - ptr);
	if (len > 2) {
	  snew(buf,(len+1));
	  strncpy(buf,ptr,len);
	  if (strstr(ptr,option[n]) != ptr) {
          set_warning_line(wi,"mdp file",-1);
          sprintf(warn_buf,"Malformed %s option %s",nopt[n],buf);
          warning(wi,warn_buf);
	  }
	  else {
	    srenew(cppopts,++ncppopts);
	    cppopts[ncppopts-1] = strdup(buf);
	  }
	  sfree(buf);
	  ptr = rptr;
	}
      }
    }
  }
  srenew(cppopts,++ncppopts);
  cppopts[ncppopts-1] = NULL;

  return cppopts;  
}


int
find_gb_bondlength(t_params *plist,int ai,int aj, real *length)
{
    int i,j,a1,a2;
    
    int found=0;
    int status;
    
    for(i=0;i<F_NRE && !found;i++)
    {
        if(IS_CHEMBOND(i))
        {
            for(j=0;j<plist[i].nr; j++)
            {
                a1=plist[i].param[j].a[0];
                a2=plist[i].param[j].a[1];
                
                if( (a1==ai && a2==aj) || (a1==aj && a2==ai))
                {
                    /* Equilibrium bond distance */
                    *length=plist[i].param[j].c[0];
                    found=1;
                }
            }
        }
    }
    status = !found;
    
    return status;
}


int
find_gb_anglelength(t_params *plist,int ai,int ak, real *length)
{
    int i,j,a1,a2,a3;
    real r12,r23,a123;
    int found=0;
    int status,status1,status2;
    
    r12 = r23 = 0;
    
    for(i=0;i<F_NRE && !found;i++)
    {
        if(IS_ANGLE(i))
        {
            for(j=0;j<plist[i].nr; j++)
            {
                a1=plist[i].param[j].a[0];
                a2=plist[i].param[j].a[1];
                a3=plist[i].param[j].a[2];

                /* We dont care what the middle atom is, but use it below */
                if( (a1==ai && a3==ak) || (a1==ak && a3==ai) )
                {                    
                    /* Equilibrium bond distance */
                    a123 = plist[i].param[j].c[0];
                    /* Use middle atom to find reference distances r12 and r23 */
                    status1 = find_gb_bondlength(plist,a1,a2,&r12);
                    status2 = find_gb_bondlength(plist,a2,a3,&r23);
                    
                    if(status1==0 && status2==0)
                    {
                        /* cosine theorem to get r13 */
                        *length=sqrt(r12*r12+r23*r23-(2*r12*r23*cos(a123/RAD2DEG)));
                        found=1;
                    }
                }
            }
        }
    }
    status = !found;
    
    return status;
}

int 
generate_gb_exclusion_interactions(t_molinfo *mi,gpp_atomtype_t atype,t_nextnb *nnb)
{
	int          i,j,k,n,ai,aj,ti,tj;
    int          n12,n13,n14;
    int          ftype;
    t_param      param;
    t_params *   plist;
    t_atoms *    at;
    real         radiusi,radiusj;
    real         gb_radiusi,gb_radiusj;
    real         param_c2,param_c4;
    real         distance;
    
    plist = mi->plist;
    at    = &mi->atoms;
    
    for(n=1;n<=nnb->nrex;n++)
    {
        switch(n)
        {
            case 1: 
                ftype=F_GB12;
                param_c2 = STILL_P2;
                param_c4 = 0.8875;
                break;
            case 2:
                ftype=F_GB13;
                param_c2 = STILL_P3;
                param_c4 = 0.3516;
                break;
            default:
                /* Put all higher-order exclusions into 1,4 list so we dont miss them */
                ftype=F_GB14;
                param_c2 = STILL_P3;
                param_c4 = 0.3516;
                break;
        }
        
        for(ai=0;ai<nnb->nr;ai++)
        {
            ti         = at->atom[ai].type;
            radiusi    = get_atomtype_radius(ti,atype);
            gb_radiusi = get_atomtype_gb_radius(ti,atype);
            
            for(j=0;j<nnb->nrexcl[ai][n];j++)
            {
                aj = nnb->a[ai][n][j];
                
                /* Only add the interactions once */
                if(aj>ai)
                {
                    tj = at->atom[aj].type;
                    radiusj    = get_atomtype_radius(tj,atype);
                    gb_radiusj = get_atomtype_gb_radius(tj,atype);
                    
                   /* There is an exclusion of type "ftype" between atoms ai and aj */
                    param.a[0] = ai;
                    param.a[1] = aj;

                    /* Reference distance, not used for 1-4 interactions */
                    switch(ftype)
                    {
                        case F_GB12:
                            if(find_gb_bondlength(plist,ai,aj,&distance)!=0)
                            {
                                gmx_fatal(FARGS,"Cannot find bond length for atoms %d-%d",ai,aj);
                            }
                            break;
                        case F_GB13:
                            if(find_gb_anglelength(plist,ai,aj,&distance)!=0)
                            {
                                gmx_fatal(FARGS,"Cannot find length for atoms %d-%d involved in angle",ai,aj);
                            }                            
                            break;
                        default:
                            distance=-1;
                            break;
                    }
                    /* Assign GB parameters */
                    /* Sum of radii */
                    param.c[0] = radiusi+radiusj;
                    /* Reference distance distance */
                    param.c[1] = distance;
                    /* Still parameter */
                    param.c[2] = param_c2;
                    /* GB radius */
                    param.c[3] = gb_radiusi+gb_radiusj;
                    /* Parameter */
                    param.c[4] = param_c4;
                    
                    /* Add it to the parameter list */
                    add_param_to_list(&plist[ftype],&param);
                }
            }
        }
	}
	return 0;
}




static char **read_topol(const char *infile,const char *outfile,
                         const char *define,const char *include,
                         t_symtab    *symtab,
                         gpp_atomtype_t atype,
                         int         *nrmols,
                         t_molinfo   **molinfo,
                         t_params    plist[],
                         int         *combination_rule,
                         double      *reppow,
                         t_gromppopts *opts,
                         real        *fudgeQQ,
                         int         *nmolblock,
                         gmx_molblock_t **molblock,
                         gmx_bool        bFEP,
                         gmx_bool        bGenborn,
                         gmx_bool        bZero,
                         gmx_bool        bVerbose,
                         warninp_t   wi)
{
  FILE       *out;
  int        i,sl,nb_funct,comb;
  char       *pline=NULL,**title=NULL;
  char       line[STRLEN],errbuf[256],comb_str[256],nb_str[256];
  char       genpairs[32];
  char       *dirstr,*dummy2;
  int        nrcopies,nmol,nmolb=0,nscan,ncombs,ncopy;
  double     fLJ,fQQ,fPOW;
  gmx_molblock_t *molb=NULL;
  t_topology *block=NULL;
  t_molinfo  *mi0=NULL;
  DirStack   *DS;
  directive  d,newd;
  t_nbparam  **nbparam,**pair;
  t_block2   *block2;
  real       fudgeLJ=-1;    /* Multiplication factor to generate 1-4 from LJ */
  gmx_bool       bReadDefaults,bReadMolType,bGenPairs,bWarn_copy_A_B;
  double     qt=0,qBt=0; /* total charge */
  t_bond_atomtype batype;
  int        lastcg=-1;
  int        dcatt=-1,nmol_couple;
  /* File handling variables */
  int        status,done;
  gmx_cpp_t  handle;
  char     *tmp_line=NULL;
  char       warn_buf[STRLEN];

  /* We need to open the output file before opening the input file,
   * because cpp_open_file can change the current working directory.
   */
  if (outfile) {
    out = gmx_fio_fopen(outfile,"w");
  } else {
    out = NULL;
  }

  /* open input file */
  status = cpp_open_file(infile,&handle,cpp_opts(define,include,infile,wi));
  if (status != 0) 
    gmx_fatal(FARGS,cpp_error(&handle,status));

  /* some local variables */
  DS_Init(&DS);			/* directive stack			 */
  nmol     = 0;			/* no molecules yet...			 */
  d        = d_invalid;         /* first thing should be a directive 	 */
  nbparam  = NULL;              /* The temporary non-bonded matrix       */
  pair     = NULL;              /* The temporary pair interaction matrix */
  block2   = NULL;		/* the extra exclusions			 */
  nb_funct = F_LJ;
  *reppow  = 12.0;              /* Default value for repulsion power     */
  
  comb     = 0;
  
  /* Init the number of CMAP torsion angles  and grid spacing */
  plist->grid_spacing = 0;
  plist->nc           = 0;
	
  bWarn_copy_A_B = bFEP;

  batype = init_bond_atomtype();
  /* parse the actual file */
  bReadDefaults = FALSE;
  bGenPairs     = FALSE;
  bReadMolType  = FALSE;
  nmol_couple = 0;
  
  do {
    status = cpp_read_line(&handle,STRLEN,line);
    done = (status == eCPP_EOF);
    if (!done) {
      if (status != eCPP_OK)
	gmx_fatal(FARGS,cpp_error(&handle,status));
      else if (out)
	fprintf(out,"%s\n",line);

      set_warning_line(wi,cpp_cur_file(&handle),cpp_cur_linenr(&handle));
    
      pline = strdup(line);
    
	/* Strip trailing '\' from pline, if it exists */
	sl = strlen(pline);
	if ((sl > 0) && (pline[sl-1] == CONTINUE)) {
		pline[sl-1] = ' ';
	}
		
		/* build one long line from several fragments - necessary for CMAP */
		while (continuing(line))
		{
			status = cpp_read_line(&handle,STRLEN,line);
            set_warning_line(wi,cpp_cur_file(&handle),cpp_cur_linenr(&handle));
			
			/* Since we depend on the '\' being present to continue to read, we copy line                                           
			 * to a tmp string, strip the '\' from that string, and cat it to pline                                                 
			 */
			tmp_line=strdup(line);
			
			sl = strlen(tmp_line);
			if ((sl > 0) && (tmp_line[sl-1] == CONTINUE)) {
				tmp_line[sl-1] = ' ';
			}
			
			done = (status == eCPP_EOF);
			if (!done) {
				if (status != eCPP_OK)
					gmx_fatal(FARGS,cpp_error(&handle,status));
				else if (out)
					fprintf(out,"%s\n",line);
			}
			
			srenew(pline,strlen(pline)+strlen(tmp_line)+1);
			strcat(pline,tmp_line);
            sfree(tmp_line);
		}
				
      /* skip trailing and leading spaces and comment text */
      strip_comment (pline);
      trim (pline);
    
      /* if there is something left... */
      if ((int)strlen(pline) > 0) {
	if (pline[0] == OPENDIR) {
	  /* A directive on this line: copy the directive 
	   * without the brackets into dirstr, then
	   * skip spaces and tabs on either side of directive
	   */
	  dirstr = strdup((pline+1));
	  if ((dummy2 = strchr (dirstr,CLOSEDIR)) != NULL)
	    (*dummy2) = 0;
	  trim (dirstr);
	  
	  if ((newd = str2dir(dirstr)) == d_invalid) {
	    sprintf(errbuf,"Invalid directive %s",dirstr);
	    warning_error(wi,errbuf);
	  }
	  else {
	    /* Directive found */
	    if (debug) 
	      fprintf(debug,"found directive '%s'\n",dir2str(newd));
	    if (DS_Check_Order (DS,newd)) {
	      DS_Push (&DS,newd);
		d = newd;
	    }
	    else {
	      /* we should print here which directives should have
		 been present, and which actually are */
	      gmx_fatal(FARGS,"%s\nInvalid order for directive %s",
			cpp_error(&handle,eCPP_SYNTAX),dir2str(newd));
	      /* d = d_invalid; */
	    }
	  }
	  sfree(dirstr);
	}
	else if (d != d_invalid) {
	  /* Not a directive, just a plain string 
	   * use a gigantic switch to decode,
	   * if there is a valid directive!
	   */
	  switch (d) {
	  case d_defaults:
	    if (bReadDefaults)
	      gmx_fatal(FARGS,"%s\nFound a second defaults directive.\n",
			cpp_error(&handle,eCPP_SYNTAX));
	    bReadDefaults = TRUE;
	    nscan = sscanf(pline,"%s%s%s%lf%lf%lf",
			   nb_str,comb_str,genpairs,&fLJ,&fQQ,&fPOW);
	    if (nscan < 2)
            too_few(wi);
	    else {
	      bGenPairs = FALSE;
	      fudgeLJ   = 1.0;
	      *fudgeQQ  = 1.0;
	      
	      get_nbparm(nb_str,comb_str,&nb_funct,&comb,wi);
	      *combination_rule = comb;
	      if (nscan >= 3) {
		bGenPairs = (gmx_strncasecmp(genpairs,"Y",1) == 0);
		if (nb_funct != eNBF_LJ && bGenPairs) {
		  gmx_fatal(FARGS,"Generating pair parameters is only supported with LJ non-bonded interactions");
		}
	      }
	      if (nscan >= 4)
		fudgeLJ   = fLJ;
	      if (nscan >= 5)
		*fudgeQQ  = fQQ;
	      if (nscan >= 6)
		*reppow   = fPOW;
	    }
	    nb_funct = ifunc_index(d_nonbond_params,nb_funct);
	    
	    break;
	  case d_atomtypes:
	    push_at(symtab,atype,batype,pline,nb_funct,
                &nbparam,bGenPairs ? &pair : NULL,wi);
	    break;
	    
	  case d_bondtypes:
          push_bt(d,plist,2,NULL,batype,pline,wi);
	    break;
	  case d_constrainttypes:
          push_bt(d,plist,2,NULL,batype,pline,wi);
	    break;
	  case d_pairtypes:
	    if (bGenPairs)
            push_nbt(d,pair,atype,pline,F_LJ14,wi);
	    else
            push_bt(d,plist,2,atype,NULL,pline,wi);
	      break;
	  case d_angletypes:
          push_bt(d,plist,3,NULL,batype,pline,wi);
	    break;
	  case d_dihedraltypes:
	    /* Special routine that can read both 2 and 4 atom dihedral definitions. */
          push_dihedraltype(d,plist,batype,pline,wi);
	    break;

	  case d_nonbond_params:
          push_nbt(d,nbparam,atype,pline,nb_funct,wi);
	    break;
	    /*
	      case d_blocktype:
	      nblock++;
	      srenew(block,nblock);
	      srenew(blockinfo,nblock);
	      blk0=&(block[nblock-1]);
	      bi0=&(blockinfo[nblock-1]);
	      init_top(blk0);
	      init_molinfo(bi0);
	      push_molt(symtab,bi0,pline);
	      break;
	    */

	  case d_implicit_genborn_params:
          push_gb_params(atype,pline,wi);
	    break;

	  case d_implicit_surface_params:
	    gmx_fatal(FARGS,"Implicit surface directive not supported yet.");
	    break;

	  case d_cmaptypes:
          push_cmaptype(d, plist, 5, atype, batype,pline,wi);
		break;
			  
	  case d_moleculetype: {
	    if (!bReadMolType) {
	      int ntype;
          if (opts->couple_moltype != NULL &&
              (opts->couple_lam0 == ecouplamNONE ||
               opts->couple_lam0 == ecouplamQ ||
               opts->couple_lam1 == ecouplamNONE ||
               opts->couple_lam1 == ecouplamQ))
          {
              dcatt = add_atomtype_decoupled(symtab,atype,
                                             &nbparam,bGenPairs?&pair:NULL);
          }
	      ntype = get_atomtype_ntypes(atype);
	      ncombs = (ntype*(ntype+1))/2;
	      generate_nbparams(comb,nb_funct,&(plist[nb_funct]),atype,wi);
	      ncopy = copy_nbparams(nbparam,nb_funct,&(plist[nb_funct]),
				    ntype);
	      fprintf(stderr,"Generated %d of the %d non-bonded parameter combinations\n",ncombs-ncopy,ncombs);
	      free_nbparam(nbparam,ntype);
	      if (bGenPairs) {
		gen_pairs(&(plist[nb_funct]),&(plist[F_LJ14]),fudgeLJ,comb,bVerbose);
		ncopy = copy_nbparams(pair,nb_funct,&(plist[F_LJ14]),
				      ntype);
		fprintf(stderr,"Generated %d of the %d 1-4 parameter combinations\n",ncombs-ncopy,ncombs);
		free_nbparam(pair,ntype);
	      }
	      /* Copy GBSA parameters to atomtype array? */
	      
	      bReadMolType = TRUE;
	    }
	    
	    push_molt(symtab,&nmol,molinfo,pline,wi);
	    srenew(block2,nmol);
	    block2[nmol-1].nr=0;
	    mi0=&((*molinfo)[nmol-1]);
	    break;
	  }
	  case d_atoms: 
          push_atom(symtab,&(mi0->cgs),&(mi0->atoms),atype,pline,&lastcg,wi);
	    break;
	    
	  case d_pairs: 
          push_bond(d,plist,mi0->plist,&(mi0->atoms),atype,pline,FALSE,
                    bGenPairs,*fudgeQQ,bZero,&bWarn_copy_A_B,wi);
	    break;
	    
	  case d_vsites2:
	  case d_vsites3:
	  case d_vsites4:
	  case d_bonds:
	  case d_angles:
	  case d_constraints:
	  case d_settles:
	  case d_position_restraints:
	  case d_angle_restraints:
	  case d_angle_restraints_z:
	  case d_distance_restraints: 
	  case d_orientation_restraints:
	  case d_dihedral_restraints:
	  case d_dihedrals:
	  case d_polarization:
	  case d_water_polarization:
	  case d_thole_polarization:
	    push_bond(d,plist,mi0->plist,&(mi0->atoms),atype,pline,TRUE,
                  bGenPairs,*fudgeQQ,bZero,&bWarn_copy_A_B,wi);
	    break;
	  case d_cmap:
          push_cmap(d,plist,mi0->plist,&(mi0->atoms),atype,pline,
                    &bWarn_copy_A_B,wi);
		break;
			  
	  case d_vsitesn:
          push_vsitesn(d,plist,mi0->plist,&(mi0->atoms),atype,pline,wi);
	    break;
	  case d_exclusions:
	    if (!block2[nmol-1].nr)
	      init_block2(&(block2[nmol-1]),mi0->atoms.nr);
	    push_excl(pline,&(block2[nmol-1]));
	    break;
	  case d_system: 
	    trim(pline);
	    title=put_symtab(symtab,pline);
	    break;
	  case d_molecules: {
	    int  whichmol;
	    gmx_bool bCouple;
	    
	    push_mol(nmol,*molinfo,pline,&whichmol,&nrcopies,wi);
	    mi0=&((*molinfo)[whichmol]);
	    srenew(molb,nmolb+1);
	    molb[nmolb].type = whichmol;
	    molb[nmolb].nmol = nrcopies;
	    nmolb++;
	    
        bCouple = (opts->couple_moltype != NULL &&
                   (gmx_strcasecmp("system"    ,opts->couple_moltype) == 0 ||
                    gmx_strcasecmp(*(mi0->name),opts->couple_moltype) == 0));
        if (bCouple) {
            nmol_couple += nrcopies;
        }

	    if (mi0->atoms.nr == 0) {
	      gmx_fatal(FARGS,"Molecule type '%s' contains no atoms",
			*mi0->name);
	    }
	    fprintf(stderr,
		    "Excluding %d bonded neighbours molecule type '%s'\n",
		    mi0->nrexcl,*mi0->name);
	    sum_q(&mi0->atoms,nrcopies,&qt,&qBt);
	    if (!mi0->bProcessed) 
	      {
            t_nextnb nnb;
            generate_excl(mi0->nrexcl,
                          mi0->atoms.nr,
                          mi0->plist,
                          &nnb,
                          &(mi0->excls));
            merge_excl(&(mi0->excls),&(block2[whichmol]));
            done_block2(&(block2[whichmol]));
            make_shake(mi0->plist,&mi0->atoms,atype,opts->nshake);
         
            
            
            /* nnb contains information about first,2nd,3rd bonded neighbors.
             * Use this to generate GB 1-2,1-3,1-4 interactions when necessary.
             */
            if (bGenborn==TRUE)
            {
                generate_gb_exclusion_interactions(mi0,atype,&nnb);
            }
            
            done_nnb(&nnb);
            
            if (bCouple)
            {
                convert_moltype_couple(mi0,dcatt,*fudgeQQ,
                                       opts->couple_lam0,opts->couple_lam1,
                                       opts->bCoupleIntra,
                                       nb_funct,&(plist[nb_funct]));
            }
            stupid_fill_block(&mi0->mols,mi0->atoms.nr,TRUE);
            mi0->bProcessed=TRUE;
	    }
	    break;
	  }
	  default:
	    fprintf (stderr,"case: %d\n",d);
	    invalid_case();
	  }
	}
      }
      sfree(pline);
      pline=NULL;
    }
  } while (!done);
  status = cpp_close_file(&handle);
  if (status != eCPP_OK) 
    gmx_fatal(FARGS,cpp_error(&handle,status));
  if (out)
    gmx_fio_fclose(out);

  if (opts->couple_moltype) {
    if (nmol_couple == 0) {
      gmx_fatal(FARGS,"Did not find any molecules of type '%s' for coupling",
		opts->couple_moltype);
    }
    fprintf(stderr,"Coupling %d copies of molecule type '%s'\n",
            nmol_couple,opts->couple_moltype);
  }
  
  /* this is not very clean, but fixes core dump on empty system name */
  if(!title)
    title=put_symtab(symtab,"");
  if (fabs(qt) > 1e-4) {
    sprintf(warn_buf,"System has non-zero total charge: %e\n\n",qt);
    warning_note(wi,warn_buf);
  }
  if (fabs(qBt) > 1e-4 && !gmx_within_tol(qBt,qt,1e-6)) {
    sprintf(warn_buf,"State B has non-zero total charge: %e\n\n",qBt);
    warning_note(wi,warn_buf);
  }
  DS_Done (&DS);
  for(i=0; i<nmol; i++)
    done_block2(&(block2[i]));
  free(block2);
  
  done_bond_atomtype(&batype);
  
  *nrmols=nmol;
  
  *nmolblock = nmolb;
  *molblock  = molb;
  
  return title;
}

char **do_top(gmx_bool         bVerbose,
              const char   *topfile,
              const char   *topppfile,
              t_gromppopts *opts,
              gmx_bool         bZero,
              t_symtab     *symtab,
              t_params     plist[],
              int          *combination_rule,
              double       *repulsion_power,
              real         *fudgeQQ,
              gpp_atomtype_t atype,
              int          *nrmols,
              t_molinfo    **molinfo,
              t_inputrec   *ir,
              int          *nmolblock,
              gmx_molblock_t **molblock,
              gmx_bool          bGenborn,
              warninp_t     wi)
{
    /* Tmpfile might contain a long path */
    const char *tmpfile;
    char **title;
    
    if (topppfile)
    {
        tmpfile = topppfile;
    }
    else
    {
        tmpfile = NULL;
    }
    
    if (bVerbose)
    {
        printf("processing topology...\n");
    }
    title = read_topol(topfile,tmpfile,opts->define,opts->include,
                       symtab,atype,nrmols,molinfo,
                       plist,combination_rule,repulsion_power,
                       opts,fudgeQQ,nmolblock,molblock,
                       ir->efep!=efepNO,bGenborn,bZero,bVerbose,
                       wi);
    if ((*combination_rule != eCOMB_GEOMETRIC) && 
        (ir->vdwtype == evdwUSER))
    {
        warning(wi,"Using sigma/epsilon based combination rules with"
                " user supplied potential function may produce unwanted"
                " results");
    }
    
    return title;
}


static void generate_qmexcl_moltype(gmx_moltype_t *molt,unsigned char *grpnr,
				    t_inputrec *ir)
{
  /* This routine expects molt->ilist to be of size F_NRE and ordered. */

  /* generates the exclusions between the individual QM atoms, as
   * these interactions should be handled by the QM subroutines and
   * not by the gromacs routines 
   */
  int
    i,j,l,k=0,jmax,qm_max=0,qm_nr=0,nratoms=0,link_nr=0,link_max=0;
  atom_id 
    *qm_arr=NULL,*link_arr=NULL,a1,a2,a3,a4,ftype=0;
  t_blocka
    qmexcl;
  t_block2
    qmexcl2;
  gmx_bool 
    *bQMMM,*blink,bexcl;

  /* First we search and select the QM atoms in an qm_arr array that
   * we use to create the exclusions.
   *
   * we take the possibility into account that a user has defined more
   * than one QM group:
   *
   * for that we also need to do this an ugly work-about just in case
   * the QM group contains the entire system...  
   */
  fprintf(stderr,"excluding classical QM-QM interactions...\n");

  jmax = ir->opts.ngQM;

  /* we first search for all the QM atoms and put them in an array 
   */
  for(j=0;j<jmax;j++){
    for(i=0;i<molt->atoms.nr;i++){
      if(qm_nr>=qm_max){
	qm_max += 100;
	srenew(qm_arr,qm_max);
      }
      if ((grpnr ? grpnr[i] : 0) == j){
	qm_arr[qm_nr++] = i;
      }
    }
  }  
  /* bQMMM[..] is an array containin TRUE/FALSE for atoms that are
   * QM/not QM. We first set all elements to false. Afterwards we use
   * the qm_arr to change the elements corresponding to the QM atoms
   * to TRUE.  
   */
  snew(bQMMM,molt->atoms.nr);
  for (i=0;i<molt->atoms.nr;i++)
    bQMMM[i]=FALSE;
  for (i=0;i<qm_nr;i++)
    bQMMM[qm_arr[i]]=TRUE;
  
  /* We remove all bonded interactions (i.e. bonds,
   * angles, dihedrals, 1-4's), involving the QM atoms. The way they
   * are removed is as follows: if the interaction invloves 2 atoms,
   * it is removed if both atoms are QMatoms. If it involves 3 atoms,
   * it is removed if at least two of the atoms are QM atoms, if the
   * interaction involves 4 atoms, it is removed if there are at least
   * 2 QM atoms.  Since this routine is called once before any forces
   * are computed, the top->idef.il[N].iatom[] array (see idef.h) can
   * be rewritten at this poitn without any problem. 25-9-2002 */
  
  /* first check weter we already have CONNBONDS: */
  if (molt->ilist[F_CONNBONDS].nr != 0){
    fprintf(stderr,"nr. of CONNBONDS present already: %d\n",
	    molt->ilist[F_CONNBONDS].nr/3);
    ftype = molt->ilist[F_CONNBONDS].iatoms[0];
    k = molt->ilist[F_CONNBONDS].nr;
  }
  /* now we delete all bonded interactions, except the ones describing
   * a chemical bond. These are converted to CONNBONDS
   */
  for (i=0;i<F_LJ;i++){
    if(i==F_CONNBONDS)
      continue;
    nratoms = interaction_function[i].nratoms;
    j = 0;
    while (j<molt->ilist[i].nr){
      bexcl = FALSE;
      switch (nratoms){
      case 2:
	a1 = molt->ilist[i].iatoms[j+1];
	a2 = molt->ilist[i].iatoms[j+2];
	bexcl = (bQMMM[a1] && bQMMM[a2]);
	/* a bonded beteen two QM atoms will be copied to the
  	 * CONNBONDS list, for reasons mentioned above 
	 */
	if(bexcl && i<F_ANGLES){
	  srenew(molt->ilist[F_CONNBONDS].iatoms,k+3);
	  molt->ilist[F_CONNBONDS].nr         += 3; 
	  molt->ilist[F_CONNBONDS].iatoms[k++] = ftype;
	  molt->ilist[F_CONNBONDS].iatoms[k++] = a1;
	  molt->ilist[F_CONNBONDS].iatoms[k++] = a2;
	}
	break;
      case 3:
	a1 = molt->ilist[i].iatoms[j+1];
	a2 = molt->ilist[i].iatoms[j+2];
	a3 = molt->ilist[i].iatoms[j+3];
	bexcl = ((bQMMM[a1] && bQMMM[a2]) ||
		 (bQMMM[a1] && bQMMM[a3]) ||
		 (bQMMM[a2] && bQMMM[a3]));
	break;
      case 4:
	a1 = molt->ilist[i].iatoms[j+1];
	a2 = molt->ilist[i].iatoms[j+2];
	a3 = molt->ilist[i].iatoms[j+3];
	a4 = molt->ilist[i].iatoms[j+4];
	bexcl = ((bQMMM[a1] && bQMMM[a2] && bQMMM[a3]) ||
		 (bQMMM[a1] && bQMMM[a2] && bQMMM[a4]) ||
		 (bQMMM[a1] && bQMMM[a3] && bQMMM[a4]) ||
		 (bQMMM[a2] && bQMMM[a3] && bQMMM[a4]));
	break;
      default:
	gmx_fatal(FARGS,"no such bonded interactions with %d atoms\n",nratoms);
      }
      if (bexcl){
	/* since the interaction involves QM atoms, these should be
         * removed from the MM ilist 
	 */
	molt->ilist[i].nr -= (nratoms+1);
	for (l=j;l<molt->ilist[i].nr;l++)
	  molt->ilist[i].iatoms[l] = molt->ilist[i].iatoms[l+(nratoms+1)];  
      } else {
	j += nratoms+1; /* the +1 is for the functype */
      }
    }
  }
  /* Now, we search for atoms bonded to a QM atom because we also want
   * to exclude their nonbonded interactions with the QM atoms. The
   * reason for this is that this interaction is accounted for in the
   * linkatoms interaction with the QMatoms and would be counted
   * twice.  */
 
  for(i=0;i<F_NRE;i++){
    if(IS_CHEMBOND(i)){
      j=0;
      while(j<molt->ilist[i].nr){
	a1 = molt->ilist[i].iatoms[j+1];
	a2 = molt->ilist[i].iatoms[j+2];
	if((bQMMM[a1] && !bQMMM[a2])||(!bQMMM[a1] && bQMMM[a2])){
	  if(link_nr>=link_max){
	    link_max += 10;
	    srenew(link_arr,link_max);
	  }
	  if(bQMMM[a1]){
	    link_arr[link_nr++] = a2;
	  } else {
	    link_arr[link_nr++] = a1;
	  }
	}
	j+=3;
      }
    }
  }   
  snew(blink,molt->atoms.nr);   
  for (i=0;i<molt->atoms.nr;i++)
    blink[i]=FALSE;
  for (i=0;i<link_nr;i++)
    blink[link_arr[i]]=TRUE;
  /* creating the exclusion block for the QM atoms. Each QM atom has
   * as excluded elements all the other QMatoms (and itself).
   */
  qmexcl.nr = molt->atoms.nr;
  qmexcl.nra = qm_nr*(qm_nr+link_nr)+link_nr*qm_nr; 
  snew(qmexcl.index,qmexcl.nr+1);
  snew(qmexcl.a,qmexcl.nra);
  j=0;
  for(i=0;i<qmexcl.nr;i++){
    qmexcl.index[i]=j;
    if(bQMMM[i]){
      for(k=0;k<qm_nr;k++){
	qmexcl.a[k+j] = qm_arr[k];
      }
      for(k=0;k<link_nr;k++){
	qmexcl.a[qm_nr+k+j] = link_arr[k];
      }
      j+=(qm_nr+link_nr);
    }
    if(blink[i]){
      for(k=0;k<qm_nr;k++){
	qmexcl.a[k+j] = qm_arr[k];
      }
      j+=qm_nr;
    }
  }
  qmexcl.index[qmexcl.nr]=j;
  
  /* and merging with the exclusions already present in sys.
   */

  init_block2(&qmexcl2,molt->atoms.nr);
  b_to_b2(&qmexcl, &qmexcl2);
  merge_excl(&(molt->excls),&qmexcl2);
  done_block2(&qmexcl2);

  /* Finally, we also need to get rid of the pair interactions of the
   * classical atom bonded to the boundary QM atoms with the QMatoms,
   * as this interaction is already accounted for by the QM, so also
   * here we run the risk of double counting! We proceed in a similar
   * way as we did above for the other bonded interactions: */
  for (i=F_LJ14;i<F_COUL14;i++){
    nratoms = interaction_function[i].nratoms;
    j = 0;
    while (j<molt->ilist[i].nr){
      bexcl = FALSE;
      a1 = molt->ilist[i].iatoms[j+1];
      a2 = molt->ilist[i].iatoms[j+2];
      bexcl = ((bQMMM[a1] && bQMMM[a2])||
	       (blink[a1] && bQMMM[a2])||
	       (bQMMM[a1] && blink[a2]));
      if (bexcl){
	/* since the interaction involves QM atoms, these should be
         * removed from the MM ilist 
	 */
	molt->ilist[i].nr -= (nratoms+1);
	for (k=j;k<molt->ilist[i].nr;k++)
	  molt->ilist[i].iatoms[k] = molt->ilist[i].iatoms[k+(nratoms+1)];  
      } else {
	j += nratoms+1; /* the +1 is for the functype */
      }
    }
  }  
  
  free(qm_arr);  
  free(bQMMM);
  free(link_arr);
  free(blink);
} /* generate_qmexcl */ 

void generate_qmexcl(gmx_mtop_t *sys,t_inputrec *ir)
{
  /* This routine expects molt->molt[m].ilist to be of size F_NRE and ordered.
   */

  unsigned char *grpnr;
  int mb,mol,nat_mol,i;
  gmx_molblock_t *molb;
  gmx_bool bQMMM;

  grpnr = sys->groups.grpnr[egcQMMM];

  for(mb=0; mb<sys->nmolblock; mb++) {
    molb = &sys->molblock[mb];
    nat_mol = sys->moltype[molb->type].atoms.nr;
    for(mol=0; mol<molb->nmol; mol++) {
      bQMMM = FALSE;
      for(i=0; i<nat_mol; i++) {
	if ((grpnr ? grpnr[i] : 0) < ir->opts.ngQM) {
	  bQMMM = TRUE;
	}
      }
      if (bQMMM) {
	if (molb->nmol > 1) {
	  /* We need to split this molblock */
	  if (mol > 0) {
	    /* Split the molblock at this molecule */
	    sys->nmolblock++;
	    srenew(sys->molblock,sys->nmolblock);
	    for(i=mb; i<sys->nmolblock-1; i++) {
	      sys->molblock[i+1] = sys->molblock[i];
	    }
	    sys->molblock[mb  ].nmol  = mol;
	    sys->molblock[mb+1].nmol -= mol;
	    mb++;
	    molb = &sys->molblock[mb];
	  }
	  if (molb->nmol > 1) {
	    /* Split the molblock after this molecule */
	    sys->nmolblock++;
	    srenew(sys->molblock,sys->nmolblock);
	    for(i=mb; i<sys->nmolblock-1; i++) {
	      sys->molblock[i+1] = sys->molblock[i];
	    }
	    sys->molblock[mb  ].nmol  = 1;
	    sys->molblock[mb+1].nmol -= 1;
	  }

	  /* Add a moltype for the QMMM molecule */
	  sys->nmoltype++;
	  srenew(sys->moltype,sys->nmoltype);
	  /* Copy the moltype struct */
	  sys->moltype[sys->nmoltype-1] = sys->moltype[molb->type];
	  /* Copy the exclusions to a new array, since this is the only
	   * thing that needs to be modified for QMMM.
	   */
	  copy_blocka(&sys->moltype[molb->type     ].excls,
		      &sys->moltype[sys->nmoltype-1].excls);
	  /* Set the molecule type for the QMMM molblock */
	  molb->type = sys->nmoltype - 1;
	}
	
	generate_qmexcl_moltype(&sys->moltype[molb->type],grpnr,ir);
      }
      if (grpnr) {
	grpnr += nat_mol;
      }
    }
  }
}
